function S = mp(Phi, Y, k)

% S = mp(, Y, k)
% M is a N by M sensing matrix.
% Y is an observed M-dimensional signal
% M < N
% k is k-sparse

%Initialize

  [m,n] = size(Phi);

  if m > n
    Phi = Phi';
  end
  A = Phi';

  N = length(A(:, 1));
  r = Y;

  S_hat = zeros(N,1);

  for t = 1:k
    % find the max inner product and its index between row vector of Phi and r
    ip = A * r;
    absip = abs(ip);
    [~,i] = max(absip);

    % project r to the row vector which is the max inner product, 
    % and the result is the i-th value of S_hat
    S_hat(i) = A(i, :) * r;
    
    % update
    r = Y - A' * S_hat;
  end

  S = S_hat;